# wczytanie pakietu
library("terra")
W pierwszym kroku musimy wylistować dane (rastry), które zamierzamy
wczytać. W tym celu możemy wykorzystać funkcję
list.files(), która jako argument przyjmuje ścieżkę do
folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy
wczytać (pattern = "\\.TIF$") oraz zwrócić pełne ścieżki do
plików (full.names = TRUE).
# listowanie plikóW z katalogu
files = list.files("dane/landsat", pattern = "\\.TIF$", full.names = TRUE)
files
## [1] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF"
## [2] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF"
## [3] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF"
## [4] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4.TIF"
## [5] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5.TIF"
## [6] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6.TIF"
## [7] "dane/landsat/LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7.TIF"
Kiedy utworzyliśmy już listę plikóW, możemy je wczytać przy pomocy
funkcji rast() z pakietu {terra} i następnie
wyświetlić metadane.
# wczytanie danych rastrowych
landsat = rast(files)
landsat # odwołanie się do obiektu wyświetla metadane
## class : SpatRaster
## dimensions : 1853, 2846, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 602535, 687915, 5742525, 5798115 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)
## sources : LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1.TIF
## LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2.TIF
## LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3.TIF
## ... and 4 more source(s)
## names : LC08_~SR_B1, LC08_~SR_B2, LC08_~SR_B3, LC08_~SR_B4, LC08_~SR_B5, LC08_~SR_B6, ...
## min values : 21, 32, 1321, 1924, 6827, 6982, ...
## max values : 60441, 61519, 65451, 65192, 63873, 65454, ...
Możemy również skrócić lub zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.
names(landsat) # nazwy oryginalne
## [1] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B1"
## [2] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B2"
## [3] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B3"
## [4] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B4"
## [5] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B5"
## [6] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B6"
## [7] "LC08_L2SP_190024_20200418_20200822_02_T1_SR_B7"
names(landsat) = paste0("B", 1:7) # skrócenie nazw
names(landsat) # nowe nazwy
## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"
# zamiana nazw
# names(landsat) = c("Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2")
Wczytanie danych wektorowych odbywa się w analogiczny sposób za
pomocą funkcji vect().
# wczytanie danych wektorowych
poly = vect("dane/powiat_sremski.gpkg")
poly
## class : SpatVector
## geometry : polygons
## dimensions : 1, 1 (geometries, attributes)
## extent : 622510.6, 659710.4, 5754795, 5784772 (xmin, xmax, ymin, ymax)
## source : powiat_sremski.gpkg
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633)
## names : TERYT
## type : <chr>
## values : 3026
Teraz możemy przygotować prostą wizualizację jednego kanału (np. ultra blue; B1) oraz poligonu.
# wizualizacja
plot(landsat[[1]])
plot(poly, add = TRUE)
Zasięg naszego obszaru analizy ograniczony jest do powiatu
śremskiego, natomiast scena satelitarna ma o wiele większy zasięg. W
takiej sytuacji możemy dociąć rastry, dzięki czemu przetwarzanie będzie
szybsze, a finalna mapa zajmie mniej miejsca na dysku. Do docinania
rastrów służy funkcja crop() i jako argumenty musimy podać
raster oraz wektor.
landsat = crop(landsat, poly)
plot(landsat[[1]])
plot(poly, add = TRUE)
Obszar rastra zmniejszył się. Jednak możemy zauważyć, że poza
poligonem wartości nie zostały usunięte. Wynika to z faktu, że obraz
zawsze docinany jest do obwiedni (bounding box), a wartości
poza obiektem/poligonem w rzeczywistości są maskowane (tj. są oznaczane
jako brakujące wartości). Aby zamaskować piksele poza poligonem należy
użyć funkcji mask().
landsat = mask(landsat, poly)
plot(landsat[[1]])
plot(poly, add = TRUE)
Tą operację można przeprowadzić również w jednej linii kodu używając
argumentu mask = TRUE w funkcji crop().
crop(landsat, poly, mask = TRUE)
W następnym kroku możemy w prosty sposób sprawdzić statystyki opisowe naszego zbioru danych.
summary(landsat)
## B1 B2 B3 B4
## Min. : 6577 Min. : 32 Min. : 6844 Min. : 7686
## 1st Qu.: 8039 1st Qu.: 8220 1st Qu.: 9093 1st Qu.: 8649
## Median : 8425 Median : 8614 Median : 9609 Median : 9353
## Mean : 8650 Mean : 8921 Mean : 9903 Mean : 9961
## 3rd Qu.: 9044 3rd Qu.: 9370 3rd Qu.:10396 3rd Qu.:10796
## Max. :19736 Max. :21314 Max. :24431 Max. :27674
## NA's :48391 NA's :48390 NA's :48390 NA's :48390
## B5 B6 B7
## Min. : 7240 Min. : 7378 Min. : 7357
## 1st Qu.:14669 1st Qu.:12482 1st Qu.:10002
## Median :17462 Median :14082 Median :11417
## Mean :18004 Mean :14673 Mean :12519
## 3rd Qu.:21254 3rd Qu.:16616 3rd Qu.:14170
## Max. :31081 Max. :50865 Max. :60130
## NA's :48390 NA's :48390 NA's :48390
Jak możemy zauważyć wartości odbicia spektralnego dla naszego zbioru danych są w zakresie od kilku do kilkunastu tysięcy dla każdego kanału. Z definicji odbicie spektralne jest w przedziale od 0 do 1. Takie dane musimy przeskalować za pomocą poniższego równania:
\[x = x \cdot 0.0000275 - 0.2\]
Nie ma potrzeby stosowania tego wzoru osobno dla każdego kanału w pętli, ponieważ domyślnie operacje matematyczne są stosowane dla wszystkich kanałów.
landsat = landsat * 2.75e-05 - 0.2
summary(landsat)
## B1 B2 B3 B4
## Min. :-0.02 Min. :-0.20 Min. :-0.01 Min. :0.01
## 1st Qu.: 0.02 1st Qu.: 0.03 1st Qu.: 0.05 1st Qu.:0.04
## Median : 0.03 Median : 0.04 Median : 0.06 Median :0.06
## Mean : 0.04 Mean : 0.05 Mean : 0.07 Mean :0.07
## 3rd Qu.: 0.05 3rd Qu.: 0.06 3rd Qu.: 0.09 3rd Qu.:0.10
## Max. : 0.34 Max. : 0.39 Max. : 0.47 Max. :0.56
## NA's :48391 NA's :48390 NA's :48390 NA's :48390
## B5 B6 B7
## Min. :0.00 Min. :0.00 Min. :0.00
## 1st Qu.:0.20 1st Qu.:0.14 1st Qu.:0.08
## Median :0.28 Median :0.19 Median :0.11
## Mean :0.30 Mean :0.20 Mean :0.14
## 3rd Qu.:0.38 3rd Qu.:0.26 3rd Qu.:0.19
## Max. :0.65 Max. :1.20 Max. :1.45
## NA's :48390 NA's :48390 NA's :48390
Nadal możemy zauważyć, że pewne wartości przekraczają nasz zakres od 0 do 1. Są to wartości odstające, które zazwyczaj związane są z błędnym pomiarem lub nadmierną saturacją. Można ten problem rozwiązać na dwa sposoby:
Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.
# sposób nr 1
landsat[landsat < 0] = NA
landsat[landsat > 1] = NA
# sposób nr 2
landsat[landsat < 0] = 0
landsat[landsat > 1] = 1
Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym
przypadku zamiast funkcji plot() należy użyć funkcji
plotRGB() oraz zdefiniować kolejność kanałóW czerwonego,
zielonego oraz niebieskiego. Oprócz tego, należy podać maksymalną
wartość odbicia dla kanałów (w naszym przypadku scale = 1).
Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto
zastosować rozciągnięcie histogramów używając argumentu
stretch = "lin" lub stretch = "hist".
plotRGB(landsat, r = 4, g = 3, b = 2, scale = 1, stretch = "lin")
# wczytanie pakietu do klasteryzacji
library("cluster")
Modele do klasteryzacji (oraz klasyfikacji nadzorowanej) wymagają
macierzy lub ramki danych. Zatem dane do treningu musimy przygotować w
odpowiedni sposób. Raster można sprowadzić do macierzy przy użyciu
funkcji values().
mat = values(landsat)
nrow(mat) # wyświetla liczbę wierszy
## [1] 1238760
Za pomocą funkcji View() możemy sprawdzić jak wygląda
nasza macierz. Jak łatwo można zauważyć, mnóstwo jej wartości
oznaczonych jest jako brak danych (głównie są to wartości poza obszarem
analizy). Zazwyczaj modele nie obsługują NA, więc musimy je usunąć.
Służy do tego dedykowana funkcja na.omit().
mat_omit = na.omit(mat)
nrow(mat_omit)
## [1] 640802
Teraz przejdziemy do najważniejszego etapu analizy, czyli do
wytrenowania modelu. Użyjemy najprostszego modelu grupowania
k-średnich (k-means). Ten model wymaga jedynie, aby podać z
góry liczbę klastrów/grup (argument centers). Jest to
algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby
analiza była powtarzalna musimy ustawić ziarno losowości –
set.seed().
set.seed(1)
mdl = kmeans(mat_omit, centers = 5)
W wyniku powyższej operacji otrzymaliśmy m. in.:
mdl$centers).mdl$cluster).Wyświetlmy te obiekty:
mdl$centers
## B1 B2 B3 B4 B5 B6 B7
## 1 0.02308493 0.02740529 0.04435163 0.04147431 0.1727330 0.1309803 0.07958614
## 2 0.02044635 0.02579641 0.05547638 0.03750665 0.4421821 0.1492823 0.07336897
## 3 0.08051378 0.09522015 0.13257100 0.16603768 0.2716112 0.3437177 0.32147881
## 4 0.03431384 0.04119050 0.07134604 0.06554308 0.3509450 0.2042120 0.12686408
## 5 0.04921918 0.05830302 0.08442430 0.09825275 0.2417769 0.2551757 0.19605182
head(mdl$cluster) # wyświetla pierwsze 6 elementów
## [1] 4 5 5 1 1 1
Oznacza to, że pierwszy wiersz (piksel) należy do grupy 3, drugi do grupy 2, trzeci do grupy 2, itd. Kolejnym etapem jest stworzenie mapy na podstawie otrzymanego wektora z klastrami.
Na początku musimy przygotować pusty wektor składający się z
całkowitej liczby pikseli rastra. Można to sprawdzić za pomocą funkcji
ncell(). W naszym przypadku jest to 1 238 760.
vec = rep(NA, ncell(landsat)) # przygotuj pusty wektor
Następnie musimy przypisać nasze klastry w wektorze w odpowiednie
miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych
wartości można odwołać się przez funkcję
complete.cases().
# zastąp tylko te wartości, które nie są NA
vec[complete.cases(mat)] = mdl$cluster
W ostatnim kroku skopiujemy metadane obiektu landsat,
ale tylko z jedną warstwą i przypiszemy mu wartości wektora
vec.
clustering = rast(landsat, nlyrs = 1, vals = vec)
Sprawdźmy teraz jak wygląda klasteryzacja naszego obszaru na mapie.
colors = c("#086209", "#fdd327", "#d9d9d9", "#29a329", "#91632b")
category = c("lasy/woda", "pola uprawne", "odkryta gleba", "roślinność", "nieużytki")
plot(clustering, col = colors, type = "classes", levels = category)
Jeśli wynik jest satysfakcjonujący, to możemy go zapisać używając
funkcji writeRaster() i wczytać w innym programie.
writeRaster(clustering, "clustering.tif")
Do analizy właściwości klastrów, zamiast statystyk opisowych, mogą zostać wykorzystane wizualizacje. Największe możliwości dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”.
ggplot2 wymaga przygotowania zbioru danych do odpowiedniej postaci. Dane muszą być przedstawione w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej zmiany można dokonać w prosty sposób używając pakietu tidyr.
# install.packages(c("tidyr", "ggplot2"))
library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych
Nasz zbiór danych jest całkiem pokaźny (blisko 9 mln wartości). Nie
ma potrzeby przedstawiania wszystkich danych na wykresie. Wymaga to
więcej RAM i znacząco wydłuża czas rysowania. Prawie identyczny efekt
można uzyskać wykorzystując mniejszą próbkę danych. Jako przykład
zobrazujmy jedynie 10 000 wartości z każdego kanału spektralnego. Do
stworzenia losowej próby służy funkcja sample(). W wyniku
otrzymamy indeksy wylosowanych wierszy.
idx = sample(1:nrow(mat_omit), size = 10000)
head(idx) # wyświetl 6 pierwszy indeksów
## [1] 294762 392686 640775 538191 270373 549593
Połączmy teraz wylosowane wiersze z macierzy z odpowiednimi klastrami
(cbind()). Następnie macierz zamienimy na ramkę danych
(as.data.frame()).
stats = cbind(mat_omit[idx, ], cluster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)
## B1 B2 B3 B4 B5 B6 B7 cluster
## 1 0.0215950 0.0250875 0.0423575 0.0363350 0.1958350 0.1312650 0.0719200 1
## 2 0.0209625 0.0296525 0.0567400 0.0480775 0.3601475 0.1481500 0.0784650 4
## 3 0.0611400 0.0707100 0.1034350 0.1209800 0.3034425 0.3146350 0.2814150 3
## 4 0.0356200 0.0396350 0.0560250 0.0620475 0.1954225 0.2082925 0.1317050 1
## 5 0.0338875 0.0392775 0.0642200 0.0594900 0.3336650 0.2240500 0.1200450 4
## 6 0.0487650 0.0554750 0.0807475 0.0867425 0.2840550 0.2358200 0.1891525 5
Jak można zauważyć, powyższe dane mają formę szeroką (każdy kanał
spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić
formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu
wykorzystamy funkcję pivot_longer().
stats = pivot_longer(stats, cols = 1:7, names_to = "band", values_to = "value")
Dla formalności możemy jeszcze zmienić typ liczbowy danych (klastrów i kanałów) na kategoryczny.
stats$cluster = as.factor(stats$cluster)
stats$band = as.factor(stats$band)
Struktura danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.
ggplot(stats, aes(x = band, y = value, color = cluster)) +
geom_boxplot()
Zmieńmy kilka domyślnych parametrów żeby poprawić stylistykę ryciny.
ggplot(stats, aes(x = band, y = value, color = cluster)) +
geom_boxplot(show.legend = FALSE) +
scale_colour_manual(values = colors) +
facet_wrap(vars(cluster)) +
xlab("Kanał") +
ylab("Odbicie") +
theme_light()
Zmieniając facet_wrap(vars(cluster)) na
facet_wrap(vars(band)), zamiast zestawienia kanałów w
poszczególnych panelach, możemy zestawić klastry.